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The Stieltjes Constants are the coefficients in the Laurent expansion of the Riemann Zeta function 
g(z) about its simple pole at z=l. They can be represented as the limit of the difference between the 
sum of the first n terms of a series and the integral of its nth. term. 

The first 20 coefficients have been computed to 15 D using the Euler-Maclaurin method. As a 
by-product the sums of the series 

r,-y (-u*oog*)"/* 

have been obtained to 15 D for n = 1 ( 1 ) 20. 
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1 . Introduction. 

The Riemann ^-function defined by 

CM = J it- (1.0) 



is regular if x = S^z > 1; a continuation into the half plane Sftz > is given by 

(2i-*-lKU) = J(-l)»n-*- (1.0') 

i 

We can avoid the introduction of the idea of continuation by defining £(z) in Sftz > 0, except 
at z=l, by (1.0'). Uniform convergence in fftz > follows since 2( — \) n n~ x is convergent as an 
alternating series for x >0. (Cf. Knopp. [10, p. 441] 1 .) The fact that 

(z-l){(*)-M 

00 

as z^> 1 then follows since V (— \) n n~ 1 = — log 2 and since the exponential series gives 

i 
(2 1 ~ z — 1) = e^~ z)log2 — 1 

— U -l)h,,[l + l!=lJ[iSi + ...]. 

This function has a Laurenl expansion about the simple pole at z= 1: 

az)=-^-r + A + A l (z-l)+A 2 (z-l) 2 + A 3 (z-iy+. . . (1.1) 

z — 1 
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where the coefficients are given by 



A n = (-l) n y n /n\ 



where 

y n = lim 



« (log/;)* (\ogm)*+i ) 

2-1 ^T~ J^-0,1,2,. . .. (1.2) 



A-=l 



(We have changed the notation of Stieltjes and Jensen to agree with that of Hardy and Briggs 
and Chowla.) 

It is clear that yo — 7, Euler's constant. Stieltjes proposed [3, I, Letter 77] the computation of 
(the first five of) these constants and gave [3, I, Letter 71] the value 

-A x = y x = -0.07281 5520 
and estimated 

A 2 = -0.0047 

Hermite [3,1, Letter 74] wrote that, during a session of the Academy of Sciences, he found some 
objections to Stieltjes' proof of (1.2) and that he had obtained a more correct proof. Stieltjes 
[3,1, Letter 75] gave a detailed proof of (1.2). This result was also obtained by Jensen ([9] and [3, II, 
p. 451]) and stated by Hardy [8] and by Ramanujan [12]. Two proofs have been given recently by 
Briggs and Chowla [2]— we reproduce one of these in section 4. 

We note that tables of the A's have been prepared by Jensen [9], who gives Au i = 1 (1) 9 to 
9Z), and by Gram [7], who gives A^i— 1(1) 16 to 16D. In section 9 we describe Gram's method of 
computation. In each table the later coefficients are zero to the number of decimals given. We 
shall see that these results appear to be correct. It is interesting to note that one of the reasons 
for calculating these constants was the determination of the small complex zeros of £(s). (See 
also Lammel [21].) The apparently irregular behavior of these coefficients has been investigated 
analytically by Briggs [1]. 

Our first calculations were based on extensions by Hardy [8] of a result of Vacca [15] which 
leads to expressions for the y's in terms of elementary constants and the series 

- k (log/r)» 

Tn — 2j \ *-) — 7 — » /I — 0,1,2, . . .. 

A=2 K 

For instance, Hardy gave 

ro= 1 -log 2 = 1-0.69314 718 = 0.30685 282 

n = - 1 (log 2) 2 + y log 2 = 0. 15986 890 

t 2 = - |(log 2) 3 + y (log 2) 2 + 2yi log 2 = 0.06537 259. 

(We have corrected errors of sign in Hardy's expressions for n, T2.) The general form of the relations 
connecting the r's and y's has been given by Briggs and Chowla and is discussed in section 3. 
We calculated the r's by use of a delayed Euler transform and then obtained the y's by solving 
the triangular system of linear equations. However (see sec. 5) this method proved unsatisfactory 
and we made direct calculations of the y's using the Euler-Maclaurin series. This is a natural 
extension of the calculations made by Knuth [11] in which he obtained y to 1271 D. Error estimates 
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are easier to obtain in the Euler-Maclaurin computation, but these computations are more lengthy 
than the direct ones. 

We have only given error estimates for the cases of yi, y 2 but representative computations 
for the later y's suggest that our values of y are secure. We have therefore to regard the y's as 
our primary results and the t's as derived from them. 

Our interest in the present problem was aroused by an entry (giving yi) in the tables of Wheelon 
[6] which we used for exercises in courses dealing with the Euler transforms. 

Professor H. Zassenhaus has called our attention to rapidly convergent series for y obtained 
by Jacobstahl [24] (see also Addison [25]) which might lead to alternative attacks on the higher y's. 

2. The y's 

The limits by which the y's are defined clearly exist, e.g., in view of a supplement to the 
integral test (Knopp, [10, p. 295]), since 

f( . (log*)* 

is ultimately decreasing because 

, (log*)*- 1 ,, 
J k (x) = [k - log xU 

which is negative for x > e k . The sequences y^\ y^\ y^ 3) , . . . where 

K ^ v k~r 1 

V = 1 

are ultimately monotone decreasing since 

7 (n)_ y (n + l) = 1 (log (n+ 1))fc + . _^_ (log „ )k+ , _.( J < 



J n 



k+1 V 6 V " k+i V 6 ' «+l 

fk(t)dt-Mn+l) (2.1) 

> iifk{t) is decreasing for t 5* n, i.e., if n ^ e k . 
We shall now investigate the size of y ( k n \ i.e., the speed of convergence of 



v=l 



to its limit y/c. It will appear that 



Ti" ) -|^ S ^+0(»-«+«) (2.2) 



which means that direct computation of yu does not come in question. We follow the analysis 
given in the case k = by Francis and Littlewood [5, (Question B7, p. 2, 19)]. 

We require the following lemma which is a simple case of the Euler-Mackurin formula. (See 
Hardy [22, p. 300, Ex. 4]). 
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Lemma. If g is three times dijferentiable in [0, 1] then 



g (D = g (0) +| [ g '(0) + b'U)] -™ g'"(o) 



where < < 1. 
PROOF. Consider 



G(x)=g(x)-g(0) -±x[g'(0)+g'(x)] -x* 



g(l )- g(0) - | (s'(O) +S'U)} 



We have 



G'W = g'W-|k'(O) + g'W]-|v'W-^(i)-«(0)-^{g'(0)+«'(i)}] 



and 



C"(*) = -|v"(*)-6* 



s(l)-*(0)-± {*'«)) + *'(!)}]. 



We note that G(0) = G(l) = and so G'(d ) = for some ft, < ft < 1. Also G'(0) = 0. 
Hence G"(0) = for some ft < < ft. Since MOwe have 

*"'(•) 12[g(l) -g(0) -| {g'(0) +g'(l)}] 

which is the result required. 

Applying the lemma in the case 



gives (cf. (2.1)): 



g(x) - (log (* +«))*+ '/(A +1), 



(log(n+l))*+' Qogn)*" (lo g (ra+l)) fc ^ (b+1) 

n + 1 r * r * 



A+l 



*+l 



= 2 [/*(*)-/*(»+!)] + «(».*) 



where 



e(n,k) = 0((logn) k ln*). 
Summing (2.3) for n = m to m = M — 1 gives for a suitable £, 



(2.3) 



1 



y(m)-ylM) >L [fk{m) -f k (M)] ~K 2 (l0g/l)^ 3 



and 



yjW _ jiM) ^ i [/fc(m) _ /fc(M)] +K ^ (log n)fc/n 3 
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Letting M -> o° gives, for any e > 0, 



y™ = \fk(m)+0{m-^) 



which is the result (2.2) announced. 

3. Derivation of the r-y Relations 

We have seen that as n — > °° : 



» (log v) k (logft)^ 1 

2— r" = -m— h ^- h0(1) - 



Hence 



4^ (log ")* (log 2n)^ +1 

S — y — = m r-y fr + o(l). 



The binomial series for (log 2v) k = (log 2 + log v) k gives 



21 



)fc= iMi ( ' )(iog2)A "' (iog ^' 



2^ 



-1 (>- a-l^' 



(3.1) 



(3.2) 



(3.3) 



i(,*)(lo g 2)~ 



(log re)' 



+ yt+o(\) 



where we have used (3.1) in the last stage. 

Integrating (with respect to /3) the binomial expansion of (a + fl) k from to b, or otherwise, 
we find 



!(')-£ 



(q + P) k+r 



fc+1 



o ifc+1 t + l" 



We use this with 
to get 



a = log2, i = logn 



V ( k \ n ?^-* (log")' +1 _ (log2n)* +1 (log2)*+* 
Jo^/ H-l A+l A+l 



Substituting from (3.4) in (3.3) we get 



» (log2^) fc = (log2n) fc+1 _ (log2) fc+1 ^^ fk 
\ 2v A:+l £ + _ 



and then subtracting (3.5) from (3.2) gives 



_^ (~l) y - 1 (logy) fc = (log2)* +1 , m /A: 



fc + 



1 + S(J Uog2)*-'<y,+ (l), 



(3.4) 



1+2 (J (log2)^-'y,+ o(l) (3.5) 
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the last term (t= k) in the summation with respect to t cancelling the jk on the right of (3.2). We 
now let n — » °° and find 

(log2)* +1 , {£» (k 



)k+l A-l /L\ 

5— +2(j(log2)*-^. (3.6) 



For k= 1,2,.. .we find 



ri=-|(log2) 2 + ylo g 2 

72=~ f (log 2) 3 + y(log 2) 2 + 2y, log 2 

t 3 = — i(log2) 4 + y(log2) 3 + 37,(log2) 2 + 3y 2 log2. 
The result for A: = is, trivially, 

4. Proof of Representation 

If we differentiate the relation 

n=l 

A; times with respect to z we get 

g<*>(z)=(-l)*jr (-l)»(logn)*/ii* 

n=l 

and, in particular, 

gC»(l)«(-l)*r*. (4.1) 

We have 

(2i-^-l) = e -u-Dio g 2_ 1= £ ( "" 1)W( } og2) W ( 2 -l)n 



and multiplying this by 



we get 



where i4_i = l. Hence 



C(s) = (2-l)-'+£/Uz-l)» 



« faf (-l)'(log2)' . Uk 



^, a)=H ffclM,^ ,4.2, 

(=1 r - 
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We now use the representation (3.6) for t* to find from (4.1) and (4.2) 

(-«• f-^T+ | (') <l<* 2 >~r,]-« |^" -<.-,. ,4.3, 



If we use this relation for &=1, 2, 3 respectively we find 

yi=— A u 
y 2 =2lA 2 . 

These suggest the general solution 

y»= {-l) n n\A n 

which is readily verified. The term outside the summation on the left in (4.3), 

(-1)*+! Qog2)*+i/(A + D, 

is the same as the term for t = k + 1 in the summation as the right. Then the t-th term on the 
left, for t = 0, 1, . . .,4 — 1, can be identified with the (k — t)-th term on the right. 

5. Calculation of the r's by Euler Transforms 

Knowing the efficiency of the Euler transform in handling the logarithmic series, it is natural 
to apply this to calculate the r's. 

Specifically, the Euler transform of £ (— l) n /n, when there is a delay of r, is 

1 _ H _... +( _ 1 ,-,i +( _ 1K _i T [i + i._i_ + |._i^_ + .. 

Since the sequence {n' 1 } is completely monotone, all the terms in the tail of the transform are 
of one sign and we can readily estimate the error. When we take n terms in the tail this error is 

io g 2-(,(o + c(, n) )=(-i)^[^ r - r+2 , r 1 + 2 3 ;;;; ra (r+n+1) + . . .] 

and is less in absolute value than 

((r+l)2 w + 1 )- 1 . 

Various hypotheses lead to various optimal choices of r,n. But these are not critical. It has also 
been pointed out that the optimal choice of r,n did not seem to depend critically on the (alternating) 
series being summed (see, e.g. Todd, [14]). 

The Euler transform is, of course, convergent if the original series is, but improvement in 
the speed of convergence of 2 (— \) n u n has only been established in case the sequence {u n } is 

completely monotone and if (u n +ilu n ) ^ k > - (see, e.g., Knopp [10, p. 253]). 
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When u n is completely monotone the terms in the remainder are all of one sign and error 
estimates are easy to obtain. However the terms in ta are not even monotone: since 

d [(logs)* ] (log*)*-' . 1 

we expect the terms in Tk to increase up to about the [e k ] -th and then to decrease. Continuing we 
see that 

d 2 [ (log*)* ! (log*)*- 2 -»/. x 2 ,«,, ±ui m 
— - = — [2 (log x) 2 - 3k log x + k(k - l)\ 



which vanishes for 



logx = -k±- VAM-8 



so that the second difference may be expected to decrease, become negative and then increase 
to positive values. 

We have not obtained useful error bounds for the Euler transforms of the Tk (k > 1 ) . 

We shall discuss these matters in some detail elsewhere. 



6. Calculation of the y's from the t's 

In section 3 we have derived the relations: 

r, = -(log2)^V(^-M) + '§( n )(log2)"-> r ,n = l,2,3, . . ., 



(3.6) 



which we have used to compute the t's from the y's. From this triangular system of linear equations 
it is also trivial to obtain the y's, given the t's. Although the solution of a triangular system is a 
numerically stable process (see, e.g., Wilkinson [17, 247]) we have obtained the exact inverse 
of the 20X20 matrix 



^ = l a ij], a ij — 



|(;),=o,u, 

j>i 



.,i-i 



i=l,2, • • • 



so that the y's can be obtained with less rounding from the t's. 



7o 



7. /log 2 



y„_,/(log2)«- 1 



= A~ l X 



ylog2 + T,/log2 




ylo g 2 + T 2 /(log2) 2 




^ 1 log2 + T„/(log2)« 
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It is easy to find A~ l when n= 10, e.g., by inverting^ and identifying the appropriate rational 
entries inA~ l . Specifically, we find 



A* = 



1 



1 


1 


















2 


2 
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1 


1 
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2 
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3 
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1 
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4 
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1 
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12 
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1 
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42 
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7 










1 
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1 
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o 














12 




24 




12 


2 


8 






1 
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7 
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1 






o 




— 















30 




9 




15 




3 " 


2 


9 






3 




1 




7 




3 


1 


1 


— 




o 




— 















20 




2 




10 




4 " 


~ 2 


10 



This naive approach will not work to find A~ l in the 20 X 20 case, owing to the condition of A. We 
are indebted to Morris Newman for providing us with the inverse in the general case: 



[A->h = 



if i<j 



i=l,2, . 



i^J 



where the Bk$ are the Bernoulli numbers (Z?o=l, jBi= — — , B? =—, Z?3 = 0, etc.)which are readily 
available as rational fractions or as decimals. Incidentally the characteristic vectors of A can be 
expressed as Stirling numbers of the second kind. 

7. Calculation of the y's by Euler-Maclaurin 

The Euler-Maclaurin summation formula can be written as 

X/O0-fVw*=^mp)+/(9)}+t 7 |^ T {^- 1) (g)-/ < ^ 1, (p)}+^(P,<7), (7-0) 

j=p Jp z j=\ w ) • 



where the Z?a-'s are the Bernoulli numbers, Bz= t, B4 — t^:, #6 = 7^, .... We give expressions for 



30' b 42' 



the remainder: 



Rkip, 9) = {2lc \ 1)l f 9 P*k+i(t) f** +i Ht) dt, 



(7.1a) 
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where 



» 2 sin 2r7rt 
P2fc+lW_ § (2^)^' 



(7.1b) 



where B 2 k(t) is the Bernoulli polynomial. For up-to-date discussions of these results see Ostrowski 
[19]. 

We use (7.0) in the following way. We first take p = X q = n: 

*(») -t/(/)- f"fU)dt=l\f(n) +/(D] +t T^rt/w^W -/*-»> (1)} + **U, »). 



(7.2a) 



We next assume/is such that f^^(t) — » for all relevant jli and that R* : (1, °°) exists. Letting fi— ► o° 
gives 

Knifi/(; - )-f , /(«)Al-|/(l)+2^-{-/W-«(l)}+18*U,«). (7.2b) 

Lj=i Ji -I z j=i (4/.> ! 



If we now subtract (7.2a) from (7.2b) we get 

1 



k ( 2 » ! " 



lim 9(m) = % n - ^/(n)-^ t^/W -«(n) +«*(«, oo). 



(7.3) 



This is the formula which is used to calculate the y's. We have to choose n, k so that Rkin, °°) 
is appropriately small. In the cases with which we are concerned, the series 



30 flo 



S'j-STfS./*- *") 



is easily seen to be divergent using the well known (Knopp, [10, p. 237]) estimates for B2J+2IB2J 
and the form of/^ _1) (0 when/(0 = (log t) k /t. Indeed, using the expression tor f^ r) (t) given below 
(7.11), we find 



l<j+l 



B 



■lj+2 



B 2j 

(2/+l)(2j + 2) 

4tt 2 



1 x / 2j+1) (») 

(2/+l)(2j + 2) /W-D(n) 



1 



(2;+l)(2/+2) 



(2/+l)(2/+2) 



which tends to infinity with j, ft and k being fixed. 

We shall discuss the behavior of Rk{n, °°) in the case of 71, y 2 using Ostrowski's estimate for 
(7.1b). We shall begin with an examination of y using (7.1a) which is a convenient way of finding 
an estimate which we require. 

Taking f(x) =x~ l in (7.3) we have, with the remainder in the form (7.1a), 



y-i+5+" 



■ + 



\-^ n -h\?^[i] dt 
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so that 

1+J+ • • •+- = logn + y + ;p- h \™P x (t)t-*dt. 
Z n In Jn 

By the Second Mean Value Theorem, there is an 7? ^ n such that 

\™P x {t)t-*dt = n-* Vp x {t)dt. 

Jn Jn 

Now it is known that the graph of P\(t) has the following form: 




J'* 

Jn 



It follows that 5= P, (t)dt 5* - 1/8. Hence 



logn + r + ^l+|+- • . + I^l 0g „ + r + i--J- 2 . 



(7.4) 



We shall use this result presently. 

By a similar argument, starting from 



7i : 



log 2 , log 3 



f .. . + i^_i (log ^ )2 _ip^ + r >i(0 ii 

Jn 



n 2 



In 



t 2 



dt 



we obtain 



1 n \2_l , log ^ _ log 1 , log 2 log n 

T (I 0gll ). + 7k+ - s ->- r+ - 5 - +b . .+_ 



^g^ i ri-ici 



1 /, x o , . log n 1 [1 — log n 



(7.5) 



(Cf. Boas and Wrench [20].) 

We now want to discuss the error in (7.3) in the case of yi, i.e., when/(r) = fi(t) = (logt)lt. 
We have to find the derivatives of/i (t) . 
LEMMA l.//fi(x)=(logx)/xtAen./brr ^ Owe have 



(-l) r (r') 
{?(*)= J, - [logx-d r ] 



(7.6) 



where do = 1 and for r > 



dr = 1+ I_ + . . ,+L. 



PROOF. Induction. 
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The error estimate given by Ostrowski is 

I** (», «) I « y (l -^) jfg 1 ! I/ 2 *'' (») | (7.7) 

and this is valid if / 2A) (x) is monotonic in (n, oo). In our case we shall examine when/ 2A,+ 1) (x) 
is of one sign in (n, <*>). In view of Lemma 1 we have to determine a value of x for which 

log* > G?2A:+1. 

Using the estimate (7.4) we see that we require 

1 1 



logx> log (2A;+l)+y+ 7 



2(2*+l) 8(2£+l) 2 
i.e., 

x > (2k + 1) ey exp( 2(2i + 1) - 8(2A+1)2 ) ■ 

ForA; = 10 this gives 

x > 21 X 1.7811 X 1.025 = 38.34. 
A rough estimate of the error in the case A: = 10, rc = 40 is, from (7.7), 

which is certainly negligible to the precision tabulated. 

We now turn to the case of 72. The analogue of Lemma 1 is 
LEMMA 2. 7/f 2 (x) = (log x) 2 /x then for r ^ 

f ( 2 r) M-"!^- ^~b T log x+ (log x) 2 ] (7.8) 

where ao = bo = and for r ^ 

b r+ i = b r + 2/(r+l) gmnsrb r = 2[l + 'g-+. . . + y],a r+1 = a r +b r /(r+l). (7.9) 

PROOF. Induction. 

In order to apply our estimates we have to find a value of * (r) such that 

E(r 9 x) = [a r -6 r log*+ (log x) 2 ] ^0 

if x^xo(r). From the recurrence relations defining the {a r }, {b r } we conclude that b r ~ 2 log r, 
flr ~ (log r) 2 which suggests that the critical value will be O(r). Computations also suggest that we 
may take £0— 0(r). We examine this theoretically. 
From the recurrence relations (7.9) we have 

vt 11 n tp( \ ° r 2 log x 
E(r+l 9 x)-E(r,x)= — — -^- ■ 
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Summing this we find 

E(r 9 x) = (log x) 2 -2(log x)d r + § (A«/(j+ 1)). 

We shall use our estimates (7.4), (7.5) to approximate £(r, x). We have, with errors of 0(1), 

£(r,*)=(log*)*-2(log*^^ 

= (log x) 2 -2 (log x) {logr+y}+(logr) 2 + 2ylog r. 

Trying # = 2r, i.e., log * = log r + log 2, in this we get, to the same accuracy, 

E(r 9 2r) = (log r) 2 + 2(log r) (log 2) -2 (log r)*-2(log r) (log 2 + y) 
+ (logr) 2 + 2y(logr) 
= 0. 

More precise investigations can be made but this discussion seems adequate in the context 
of practical computation. 

We shall now indicate an alternative approach. Let t r be the larger root of the quadratic 

a r — brt + t 2 = Q 9 i.e., t r = - (b r +vb 2 r — 4a r ). Then any number larger than exp t r can be taken as 

Xo. From the recurrence relations (7.9) we have, if A r — b\ — 4a r , 

Ar+i-A r =4(r+l)- 2 . 

Summing this we find 






/ 2W2 A 

<__ 4 



1 X -+. 



l(r+l)(r+2) (r + 2)(r + 3) 



2tt 2 4 



3 r+1 



From (7.4) we find 



/it 2 1 \ 1/2 

fr<logr+Y+ l + ^__L_) 

<Iogr+y+ l + ^( 1 __l_) 

<logr+( y+ ^)+l- ^ +1) 

4= log r+ 1.8599. 
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This estimate is in agreement with computations. We find 

a 25 = 0.1296 X 10 2 , 6 25 = 0.7632 X 10 1 , log 25 = 0.3219 X 10\ * 25 = 0.5083 X 10 1 

Hence we can take xo(r) = 6.5 X r. 

We give here a typical exploratory calculation which we made in the case of 72: 



n = 400: 3 Bernoulli terms give - 0. 12455 57727 57476 
which with the head of series gives y 2 = - 0.96903 63192 87234 

We use the Ostrowski error estimate (7.7) which in the present case is: 



(-4) 

(-2) 



( 1- 2V 5T ~r{ a 6~ 6 6logrc+ (logrc) 2 } 



This is valid if log n exceeds the larger root of a 6 — &6 £ + £ 2 = 0, where a 6 = 203/45 and 6 6 = 49/10, 
which is ^e = 3.6712. With n = 400, logrc == 6. The term in braces {. . .} is about 11 and so the trun- 
cation error estimate is about 10~ 19 — so that the round off errors will dominate. Note that in dealing 
with the head of the series we are summing 400 terms, each involving the square (or cube) of a 
logarithm. 

We have not made any detailed error analyses for the cases of y m , for m > 2. In our final 
computations we have taken n = 400. For instance , in the case of 715 , 

3 Bernoulli terms give 0.36042 49967 32360 (+6) 

which with the head of the series gives 715 = - 0.28369 16022 44191 (-3) 

so that there is a massive cancellation. Similarly, in the case of 719, 

3 Bernoulli terms give 0.67068 88575 83315 (+9) 

which with the head of the series gives 7i 9 = 0.21630 66613 69313 (— 3) 

It is because of this that we used the Tooper multiple precision package, with quadruple precision. 
We conclude this paragraph by recording the recurrence relations which we used to calculate 
the derivatives needed in the Bernoulli terms and in the error estimates. 
Lemma 3. //f m (x) = (log x) m /x then 

^(x)=(-l) r (r!)x- r -'[(-l) m a , r -(-l) m - i a 1 , r logx+. . . + a m , r (logx) m ]. (7.10) 



where 



«0, r+l — OLo, r + 



r+1 



ai,r+i = «i, r +; 



+ 1 



ttl, 



a.2,1 



m 



(7.11) 



r+ 1 

OLm. r — 1 
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with appropriate initial conditions: 

ato,o = tti,o = . . . = «/«-!, o = 0, a m ,o— 1. 

PROOF. Induction. 

Using the same ideas as in our discussion of / (r) we can estimate 

a*,r=M J(logr) m "* 

which gives 

/#(*) - (-D r r!^-Hlog^-logr]- 

and indicates that /J/? (*) will be of one sign for x 2= xq (r) where #o (r) = 0(r). 

8. Gram's Calculations 

Gram [7] gives few details of his computations but we can indicate the lines along which they 
went. He realized that it was not attractive to extend the calculations of Jensen [9] which were 
based on the definition 



£(5)=lim 



r=l J 



and sought another approach. 
The function 



f(*)-y*(*-l)ir-«»r(yj){(*) (8.1) 



is entire and satisfies the functional equation 
If we write 



-(*)=£(*), s = y+H (8.2) 



then = )(t) is an even entire function which does not vanish at t = since = (0)=£(— J : 
- -g- 7T- 1 / 4 T Cj\ £ (y \ We can write 



log = (it) =ao-\-ait 2 + a 2 £ i + . . . 



where ao = log £ (~5~)« If we can compute the a's, we can, by exponentiating, get the coefficients 

in the Maclaurin expansion of = (t) and from this, using (8.1) and (8.2) obtain the coefficients 

in the Laurent expansion of £ (5 ) about 5=1. 

Gram had available Stieltjes' 32D values of £(5), 5 = 2(1)70, i.e., in principle, the values of 

3 139 
log = (it) for ± it = -z-(l)-y— The a's are essentially the values of the (even) derivatives of 
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log = (it) at and can be calculated by using appropriate Lagrangian expressions. For instance, we 
have (see, e.g., [23]) 

/i/2 = ^ (3/_ 2 ~25/_ 1 4-150/o^l50/ 1 -25/ 2 + 3/ 3 )-^Mo 6 /i/2+. . . 

1 259 

hTii i =-^(-5f-2 + S9f. 1 -Wo-Wi+^f 2 -5f 3 )+^^ f j^f 1 i2-. • •• 

From these expressions we find, when/-n=/n+i — log — (i(l/2 + /i)) = log^(«+ 1), 

a =~ 0.69892 210, 

a x = + 0.02310 424. 

These are to be compared with the values given by Gram 

a = -0.69892 22679 45331 4, 

a x - + 0.0231049931 15419 0. 

Gram, actually, made use of independent calculations by himself and Stieltjes of £(i) which gave 

ao directly and his calculations were more complicated and elaborate than ours, e.g., he used the 

1 29 
values of log = (t) for it = — (1) — . 

9. Tooper's Multiple Precision Package 

When we realized that we needed multiple precision calculations we were fortunate to have 
available a package produced by Professor R. F. Tooper. In this we used one storage place for the 
exponent of our numbers and four for the mantissa (actually it is possible to use up to 19). Besides 
the usual arithmetic operations included in this package we used the logarithm subroutine. The 
logarithms were calculated from the usual series 

log~ = -2 [i + hHh 5 . . .] 
1 -r x 

after a reduction in the range to(i, 1) ; a "long division" process was used. 

Our preliminary calculations were made on the IBM 360/75 at the California Institute of 
Technology and the final ones on the IBM 370/155 at the University of South Florida. 

We used (7.3), with ra = 400, to compute the y's and A's, and (3.6) to derive the r's. 
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We are indebted to Peter Weinberger and Herman P. Robinson for detailed comments on our 
tables and especially to Henry C. Thacher, Jr. for showing us his manuscript [26] which includes a 
table of y r , r=0(l)35, 44 to 28 S as well as the coefficients of z£(l + z) together with Chebyshev 
series for this function in the ranges [—1/2, 1/2] and [0, 1], His values were also obtained by the Euler- 
Maclaurin method and our values agree with his up to the last digit except as indicated: specifically 
Thacher has for y% . . . 9|509 . . . and for y 9 . . . 0|481 . . . where we have . . . 9| and 
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